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Section  1 
INTRODUCTION 


Three-dimensional  flow  effects  play  an  important  role  in  the  perfor¬ 
mance  of  axial-flow  fans  and  compressors  that  operate  at  transonic  speeds. 

The  coupling  between  transonic  and  three-dimensional  effects  limits  the  ap¬ 
plicability  of  the  two-dimensional  analysis  methods  that  have  been  in  use  for 
some  years.  Efforts  to  extend  these  analyses  to  three-dimensional  transonic 
cases  have  been  aided  greatly  by  the  development  of  computational  methods 
for  solving  comparable  problems  in  external  aerodynamics.  The  applicable 
external-flow  methods  can  be  divided  broadly  into  two  fields:  those  based 
on  the  potential- flow  approximation  and  uuse  that  start  from  the  Euler  equa¬ 
tions.  The  potential-flow  category  is  further  divided  into  the  range  of 
small  disturbances  and  the  range  where  the  full  nonlinearity  of  the  problem 
must  be  accounted  for. 


The  nonlinear  small-disturbance  potential  theory  was  developed,  in  a 
previous  AFOSR-sponsored  study  at  Calspan.  ”  That  work  consisted  essentially 


1.  Rae,  W.J.,  "Nonlinear  Small-Disturbance  Equations  for  Three-Dimensional 
Transonic  Flow  Through  a  Compressor  Blade  Row",  AFOSR-TR-76- 1082 , 

AD-A  51254  (August  1976). 

2.  Rae,  W.J.,  "Relaxation  Solutions  for  Three-Dimensional  Transonic  Flow 

Through  a  Compressor  Blade  Row,  in  the  Nonlinear  Small-Disturbance 
Approximation",  AF0SR-TR-76-1081 ,  AD-A032555  (August  1976). 

5.  Rae,  IV. J. ,  "Finite-Difference  Calculations  of  Three-Dimensional  Transonic 
Flow  Through  a  Compressor  Blade  Row,  Using  the  Smal 1 -Disturbance  Non¬ 
linear  Potential  Equation",  pp.  228-252  of  Transonic  Flow  Problems  in 
Turbomachinery,  ed.  by  T.C.  Adamson  and  M.F.  Platcer,  Hemisphere 
Publishing  Corporation,  Washington,  (1977). 

4.  Rae,  W.J.,  "Calculations  of  Three-Dimensional  Transonic  Compressor  Flow- 

fields  by  a  Relaxation  Method",  Journal  of  Energy,  1_  (1977)  284-296. 

5.  Rae,  IV. J. ,  "Computer  Program  for  Relaxation  Solutions  of  the  Nonlinear 

Small-Disturbance  Equations  for  Transonic  Flow  in  an  Axial  Compressor 
Blade  Row",  AFOSR-TR-78-0835 ,  AD-A055744  (April  1978). 
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of  an  application  of  the  line-relaxation  methods  and  Mach-number-dependent 
differencing  procedures  pioneered  by  Murman  and  Cole.^  Flow  field  calculations 
were  done,  with  the  resulting  computer  code,  for  several  blade  rows  and 
operating  conditions.  These  calculations  showed  interesting  interactions  be¬ 
tween  the  regions  of  subsonic  and  supersonic  flow  that  develop  within  the 
blade  row. 

The  principal  limitation  of  these  results  is,  of  course,  the  small- 
disturbance  assumption.  The  pressure  ratios  and  turning  angles  of  practical 
compressors  exceed  the  range  that  can  properly  be  called  a  small  perturbation 
of  the  inlet  conditions.  Thus  the  role  of  the  previous  work  is  chiefly  to 
give  qualitative  information  about  the  flow. 

The  present  research  was  undertaken  with  the  aim  of  extending  this 
earlier  work,  so  as  to  handle  more  fully  the  nonlinearity  of  the  problem. 

As  noted  above,  the  two  principal  candidates  for  achieving  this  goal  were  the 
methods  for  solving  the  full  nonlinear  potential  equation,  and  the  time-marching 
methods  used  for  solving  the  Euler  equations.  The  former  approach  has  the 
advantage  that  only  a  single  dependent  variable  needs  to  be  stored,  compared 
with  five  dependent  variables  in  the  latter  case.  However,  the  potential-flow 
approximation  is  restricted  to  isentropic  flow.  An  additional  consideration , 
of  considerable  importance  at  the  start  of  this  research,  was  that  methods 
for  treating  the  three-dimensional  full  potential  equation  were  not  yet  de¬ 
veloped.  In  contrast,  a  number  of  papers  describing  the  "fully  implicit"  time 
marching  procedure  had  been  published,  and  appeared  to  be  capable  of  yielding 
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results  in  a  relatively  straightforward  way. 

6.  Murman,  E.M.,  and  Cole,  J.D.,  "Calculation  of  Plane  Steady  Transonic 

Flows”,  AIAA  Journal  9  (19"1)  114-121. 

7.  Beam,  R.M.  and  Warming,  R.F.,  "An  Implicit  Finite-Difference  Algorithm 
for  Hyperbolic  Systems  in  Conservation-Law  Form",  J.  Comp.  Phvs.  22_ 

(1976)  87-110.  '  ~ 

S.  Steger,  J.L.,  "Implicit  Finite  Difference  Simulation  of  Flow  About  Arbi¬ 
trary  Geometries  with  Application  to  Airfoils",  AIAA  Paper  ""-665 
(June  1977). 

9.  Kutler,  P.,  Chakravarthy ,  S.R.,  and  Lombard,  C.P.,  "Supersonic  Flow  Over 
Ablated  Nose  Tips  Using  an  Unsteady  Implicit  Numerical  Procedure", 

AIAA  Paper  78-213  (January  1978). 


Accordingly,  the  time-marching  method  was  selected  for  application 
to  the  case  of  flow  through  an  isolated  compressor  blade  row.  The  adaptations 
required  include  a  coordinate  transformation  suitable  for  a  cascade  geometry, 
modifications  to  enforce  mass-flow  conservation,  boundary  conditions  upstream 
and  downstream  of  the  blade  row,  and  a  means  of  accounting  for  the  vortex 
sheets  which  trail  downstream  of  the  blades. 

The  section  below  contains  a  review  of  the  basic  equations  in  absolute 
and  relative  coordinates,  including  two  versions  of  the  energy  equation.  Also 
described  in  this  section  are  means  for  including  the  radius  terms  that  appear 
in  cylindrical  coordinates,  and  some  details  about  the  coordinate  transforma¬ 
tion  used  and  the  metrics  that  result.  The  third  section  is  a  review  of 
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the  finite-difference  method,  patterned  after  the  Beam-Warming  technique,' 
while  the  fourth  section  contains  a  description  of  the  boundary,  wake,  Kutta, 
and  exit  conditions. 

All  of  these  elements  were  incorporated  into  a  computer  code,  and 
a  number  of  attempts  were  made  to  carry  out  a  sample  calculation.  These  ef¬ 
forts  were  not  successful,  due  principally  to  the  destabilizing  effects  of 
singularities  in  the  metric  coefficients  at  the  blade  trailing  edge,  and  at 
the  points  corresponding  to  upstream  and  downstream  infinity.  Problems  arising 
from  the  points  at  infinity  were  overcome  successfully,  but  no  satisfactory 
resolution  of  the  trailing-edge  problem  was  found. 

The  section  on  Concluding  Remarks  presents  some  suggestions  for 
further  modifications  that  may  be  capable  of  treating  the  trailing-edge  region 
successfully,  and  a  review  of  other  computer-program  elements  that  will  need 
further  development,  once  the  metric-induced  instabilities  are  removed. 
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Section  2 
BASIC  EQUATIONS 


Conservation-Law  Forms 


The  Euler  equations  for  unsteady  three-dimensional  flow  in  cylindri¬ 
cal  coordinates  x  ,  r  ,  d  may  be  written  in  conservation  form  as  (see,  for 
example.  Reference  10) 

i<7»  *  jzC-PV  *  ■fcf'-pV  »  -  * 

-~(^PV  *  j^(n  <■  yrin>vM  *  -  o 

y(rpVr)  f  J-UpVx</r )  *  ±  (r[p  +  pvr!])  +  y(PW  -  P  r  pv* 

Tt(rpV°)  +  Jx  (  rp  v*  v* )  *  7F(r°yr^)  *  je(r*Pve>  -  -  PW. 

h(re>  +  7x(r^[e +  £:(rVrC*'ri)  *  -  » 

(2-1) 


where  V  is  the  absolute  velocity,  p  and  p 
is  the  total  energy  per  unit  volume: 


are  the  pressure  and  density,  and  t 


«  *  P[CrT  +  i  l/‘] 


10,  Vinokur,  M. ,  "Conservation  Equations  of  Gasdynamics  in  Curvilinear  Co¬ 
ordinate  Svstems",  Journal  of  Computational  Physics,  _14_  (1974) 

pp.  105-125. 
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They  may  be 


These  equations  are  written  in  the  absolute  coordinates.11 
cast  in  terms  of  blade-fixed  coordinates  by  the  transformation 
t  =  t  ,  r'  =  r  ,  Q  =  6  -  cot  ,  x  '  =  X 

=  Vx  ,  Vr  >  ^  ~ 


(2-3) 


where  W  is  the  velocity  relative  to  the  blades.  After  dropping  the  primes, 
these  equations  have  the  form: 

3(pr)  %  ,  d 

Tt  +  dx  (rPVJ*'>  +  +  jQ(P*a>  -0 


at  Cr^x)  +  (rL+  +  P»*l)  +  -k  +  -L  (p^^0)  = 


3  . 


dt  Crp^r)  +  J7(rP  +  +PUrl)  +  =  P+/0(*9+u>r) 


~irp^d)  +  —(rpU%ug)  +  —(rpnruj9)  +  A  (^+p^)  = 

-1(^1,  ,  ±(r^xpl)  ,  !(«,,!)  *  (:.4) 


where  I  is  the  rothalpv: 


1  -  ST  +  2 


1  \  2 

2  (wr) 


(2-3) 


There  are  two  features  of  these  equations  that  make  it  difficult  to  apply  the 
time-marching  algorithms  developed  for  external  flow.  The  first  is  the 


11.  K’u ,  C.H.,  "A  General  Theory  of  Three-Dimensional  Flow  in  Subsonic  and 
Supersonic  Turbomachines  of  Axial  -,  Radial  -,  and  Mixed-Flow  Types”, 
N’ACA  TNT  2604  (January  1952). 
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appearance  of  the  time  derivative  on  the  right  side  of  the  energy  equation; 
in  seeking  a  steady-state  solution,  it  is  not  clear  whether  this  term  should 
be  set  equal  to  ;ero  or  evaluated  from  previous  time  steps,  and  the  external- 
flow  literature  offers  no  guidance  on  this  question.  A  preferable  form  of 
the  energy  equation  is 


Q  (rK)  +  4z  (r^CK+pl)  +  (rWr[K  +  1 p])+  (VJ&  [K  +p J  ) 


it 


2% 


=  0 


(2-6) 


where 


k  =  pi  -f  -  -rr;  *  t  [w2-  iurt!] 


(2-7) 


A  second  awkward  feature  of  these  equations  is  the  appearance  of 
the  variable  f  ,  inserted  in  several  places  in  order  to  preserve  strict  con¬ 
servation-law  form.  These  appearances  require  frequent  numerical  evaluations, 
most  of  which  can  be  avoided  if  the  strict  conservation  form  is  relaxed 
slightly,  by  associating  the  r-  factors  only  with  the  d-  derivatives: 


£  *  -kw  *  $<^*71 hw  - 


~CoVJx)  +•  J-C/p  +  pi«j2)  +  -L.(pVJxUr)+  -L  j-  (puxu&) 


~k  (/0^)+  +  ~  (f>*pU?)  +  ~ 

z 


-p  +  p  k/r  49  +  p[Wa  +  cor] 

~  “  — - -  +  — - - 


at  +  +  1-  (purw8)  +  -L  JL(p+p vi* 


dr 


P  p  Wr(atjr  +  Wg,) 


ir  +  77  (u%CK*e))  +  ^.(wrOf**o)  +  y  (w9«^.  -  !£<*♦*.) 


r  3® 


6 


(2-8) 


Dimensionless  Forms 


The  Euler  equations  can  be  made  dimensionless  in  terms  of  reference 

values  of  length,  Lref  velocity,  Uref.  and  density  prt^  ,  and  by  dividing 

2 

the  continuity,  momentum,  and  energy  equations  by  Pre^  U  ref  ,  Pref 

and  Pre.$  U  ,  respectively.  Thus  each  term  in  the  above  equations  may 

be  considered  as  a  dimensional  quantity,  or  its  dimensionless  equivalent. 

The  time,  for  example  may  be  taken  in  physical  units,  or  in  units  of  • 

These  equations  can  be  written  in  terms  of  the  following  five-compo¬ 
nent  vectors: 


where 


dU 

dt 


3£ 


3F 

2r 


1  3  Or 


H  ~  F 


r  20 


U  = 


lp  \ 

/  />w*  \ 

1  -p  4-  p  W  £  1 

P  W  r 

1  ,  e  = 

P 

p  w  ^  we 

1 

\  K  / 

\  VJx(K+p)  I 

F  = 


(2-9) 


(2-10' 


G-  = 


(  K  +  p )  / 


0 

0 


H  =  j  jo  +  p(Ws+ujry 

i  — 10  Wr  (Sujr  +  'rJg) 


(2-11) 


If  these  equations  are  subjected  to  the  general  coordinate 
transformation 

T  -  t  ,  %  =  %  (%,r,&,  t)  ,  /[  */{(X,r}&yt),  $  *  ? 

(2-12) 

1? 

then,  by  following  Viviand's  "  derivation,  it  can  be  shown  that  the  Euler 
equations  retain  conservation  form,  i.e., 


9  ,'r  lU$t+  +  F$r]  +6$gy  g  /  r  [U  gt  +  +  Qge 

3%  \  r  <8  J  9n  ^ 


_L  (  r  £UZ>t  +  E  %TL  +  F  %  rl  +&Z>Q 
9Z,  \  r  £> 


where  the  zero  term  on  the  right  side  results  from  adding  appropriate  terms 
to  achieve  the  conservation  form,  and  whereof  is  the  Jacobian: 


9  ,  n,  ?? 

9  Cx  ,  r ,  9) 


*>x  £r  ^ 

H*  Hr  % 

^  %  %  r  ^  o 


(2-14) 


These  can  also  be  written  as 


12.  Viviand,  H. ,  "Formes  Conservatives  des  Equations  de  la  Dynamique  des 
Gaz",  La  Recherche  Aerospatiale  (1974),  No.  1,  January- February , 
pp.  65-66. 
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3  0  3  £  2  F  2  Or  G 

ar  4  a^  +  T^T  *  It  =  (2‘15) 


where 


F 

A 

Gr 


H 


[U<4  +  £?„+  FKr  *  ^r}/» 
[0l1t  +  H^«  +  fr'l^t  ~]  /  £> 
[U4e  *  £<*  ^4,  *  -V2]/-©- 

H  -  F 
r£> 


(2-16) 


Finally,  by  inserting  the  definitions  of  E,  F,  and  G,  the  vectors 
r*  ^  ^ 

E  ,  F  ,  and  G  can  be  written  as 


P 


l aJ  , 


£.1 

JS> 


,0  W,  Wx  +  p 
p  W 1  Wr  +  p  ^  r 
p  W,  We  +  pEy^/r 

(K  +  p)W,-p$t  I 


^*TS 


f 

Pwz^x  +  PPx 
p  U/^  U/r  +  p  Hr 

/0  wz  +  f»  n  »/r 

(K  4  -£)  iv2  -  p  q  t 


»  oi 


/  'w*  \ 
pw3w^^t 

p  w3  U/r  -p  ZJ  r 
PW5W9  +  #>  iV'" 

(K+p)U/3-pzJ 


(2-17) 

where  Wj|  ,  and  are  the  contravariant  components  of  the  velocity  vector: 
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w,  * 

4t  + 

+ 

+ 

"e 

K  e 
r 

= 

it + 

Kcl*  + 

Ur  f\r  + 

Uq 

% 

r 

W  3  * 

+ 

+ 

5o 

r 

(2-18) 

Note  that  in  most  places  the  factor  C  appears  as  a  divisor  of  the  metrics  of 
the  third  column  of  <&,  for  example  as  £e/r;  thus  only  the  ratio  /r  needs 
to  be  stored  at  each  grid  point,  rather  than  both  factors. 
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Section  3 

FINITE-DIFFERENCE  METHOD 


The  algorithm 
(here,  the  superscript 


presented  by  Beam  and  Warming'  is  applied  as 
n  denotes  the  time  level  of  the  solution): 


follows : 


+■ 


At 

2 


H 


H 


0  (At) 


(3-1) 


Next,  a  Taylor-series  expansion  is  made,  i.e.: 


A 

-  n  + 1 
£  = 

£  n  i-  R*(Un*'- 

o") 

+  0(At)Z 

p  n  *  i 

Bn(Un"~ 

u-> 

+ 

Lr  - 

Gn  +  c"(U^- 

A  ^ 

U  ) 

+ 

D"<S"W- 

(3° ) 

+  0(At)2 

where  the  coefficients  : 

In  the  expansion  are  the 

matrices : 

•A 

A 

A 

A 

2E 

A 

£<S  A 

a  h 

R 

"  A  ) 

B  =  —  , 

c  = 

-r-sr-  ,  D  ~ 

— x- 

2  U 

2U 

3U 

at) 

A 

A  * 

A  /•  /V 

A 

A  A  A 

A-  A 

£ 

=  R  U  , 

F  -  8  0  , 

£  = 

CU  ,  H  • 

DU 

This  enables  the  equation  to  be  written  as 
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.  *£  h(>l  .  if.  +  il/  ♦  4 

2  1  Wf  *0  3<>'  H 

+  i-s"(u"'- G”;  *  Xc-cu-'-G";} 
?0  <K  J 

f  n  m  ,  *  n*i  ^  ,  1  At  ,3 

t  j  2W  +  D  (  ll  +  U  j  p  —  t  0(At) 


f  pent; 


(3-4) 


which  can  be  rearranged  as 


2  ^ 


-  At 


dn 

±Rn  + 

J-8° 

2/7 

2^ 

C-D1Y  U 


(3-5) 


where  I  is  the  identity  matrix.  Certain  terms  are  now  added,  of  order  (At) 
and  (At)3  ,  which  have  the  effect  of  "completing  the  cube"  on  the  left-hand 
side,  so  that  it  can  be  factored  as: 


-{(i*  f  f  [&c"-<])}5, 

+  o(At3) 

1  ^  ^  j 


+  —  -  H  i  +  0(At3) 


(3-6) 


Following  Beam  and  Warming,  this  equation  is  rearranged  as 
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(1  +  R'')(I+Ar-LBn)(U  ^L[4-  crt-  5” 

v  2  3£  A  2  3 n  A  2  L  it, 


AU 


_  fa£  pf  3&  r. 

=  -  AT  -s  —  +  - —  +  — -  -  H 

\3$  *n  a 


O(At)' 


(3-7) 


where 


au  =  u r,+1-  u 


The  matrices  A,  B,  C,  and  D  (without  the  (  ^  )  symbol)  are  defined  as 


ft  =  ££-  , 

s  -  i£_. 

C  «  —  , 

D  = 

3H 

au 

3  U  9 

du 

£  =  ft  U  , 

F  =  5  U  , 

Cr  -  C  L)  y 

H  = 

DU 

(3-8) 

The  relations  between 

the  two  sets  are 

A 

ft 

=  I  C  £•  + 

+  s?r  *  C%- 

A 

B 

=  1*7*  + 

*s1t  +  C  £ 

A 

c 

*  Iftt  + 

mr  »  C7 

A 

D 

D -3 

(3-9) 

r 

The  matrices  A,  B,  C, 

and  D  are  given  : 

in  the  Appendix. 
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Equation  (3-7)  is  usually  referred  to  as  the  "delta  form”  of  the 
algorithm.  Its  numerical  solution  is  found  by  a  sequence  of  three  one¬ 
dimensional  solutions: 


(J  + 

(** 
(I  + 


Here  the  term  D  has  been  placed  in  the  second  step,  in  order  to  facilitate 
the  calculation  of  two-dimensional  cases,  which  bypass  the  radial  solution 
step  altogether.  This  term  can  be  placed  in  any  one  of  the  three  steps,  with 
no  change  in  the  truncation  error. 


Damping  Terms 


The  numerical  algorithm  described  above,  which  uses  central  dif¬ 
ferences  for  the  spatial  derivatives,  requires  the  addition  of  certain  damping 
terms  for  stability.  These  terms  are  added  by  rewriting  Eq.  5-10  as  follows: 
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The  fourth  derivatives  on  the  right  hand  side  are  evaluated  explicitly,  using 
central  differences  of  data  at  the  previous  time  step.  At  grid  points  next 
to  the  boundaries,  the  second  derivative  is  used. 


The  additional  terms  on  the  left  side  are  treated  implicitly,  i.e., 
they  appear  as  corrections  to  the  coefficients  of  the  block  tridiagonal 
matrix  equations. 


Step-Size  Considerations 


For  a  fully  implicit  method,  the  time  step  would  be  limited  only  by 
considerations  of  accuracy,  and  not  by  stability.  As  will  be  noted  below, 
the  present  application  (and  all  of  the  external  -  aerodynamics  literature 
as  well)  uses  boundary  conditions  that  are  explicit,  and  involves  data  one 
time  step  behind.  This  introduces  the  problem  of  stability  considerations,  and 
it  is  usually  found  that  the  time  step  must  be  on  the  order  of  that  given  by 
the  Courant-Friedrichs-Lewy  condition: 


At 


mO-X. 


vru-Kt  c  A  §■  ,  A  q.  ,  A  <  ) 

^  maX 


(3-12) 


where  A  denotes  an  eigenvalue  of  the  matrices  that  appear  in  the  nonconserva¬ 
tive  form  of  the  Euler  equations.  For  the  present  case,  these  differ  only 
slightly  from  the  rectangular-coordinate  set  discussed  by  Warming,  Beam,  and 
Hvett.  *  The  non-conservative  forms  of  the  equations  are 


d  u. 
Tt 


+ 


d  U- 
~  ~d& 


(3-15) 


where 


13.  Warming,  R.F.,  Beam,  R.M.  and  Hyett,  B.F.,  "Diagonalization  and  Simulta¬ 
neous  Symmetrization  of  the  Gas-Dynamic  Matrices",  Mathematics  of 
Computation,  29 (1975)  1037-1045. 
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The  eigenvalues  of 

P\r  ,  and  P 

are  1 

to  those  of 

the  conservation-law  form  by 

-7  - 

P-  ^ 

where  R.  is 

Rx,  R,  ,  or  R  ,  P, 

is  R  , 

Jacobian  matrix  do/^u. 

» 

/  > 

a  o 
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i 

1  s 
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O 

P-  = 

W  y 

O  P 

O 

I 

O  0 

P 

! 

vJz-c^>‘iyi' 

\  , 

r 

c2  = 


2*> 

T- 


(3-14) 


±  c  ,  and  are  related 


P- 


(3-15) 

or  C  (see  Eq.  5-8) ,  and  u.  is  the 
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Thus 


H  A  +  A  Wx  *  ^zWr  +  U<9  j  £<>  +  +  £ZK  +  A  ^5 


±  c  V^T+^T+t (T 


(3-17) 


here  &0  ,  >  &z  >  and  $3are  either  ^  ,  % t  ,  <fr  ,  -p-  ,  or  f\t  ,  ,  f[r  , 


r  ’ 


or  ’  £<.  ’  ^r 


Coordinate  Transformations  and  Metrics 


The  geometry  of  the  axisymmetric  flow  passage  is  used  first,  to 
define  the  coordinate  i]  as  the  fractional  distance  from  hub  to  tip: 


(The  bullet-nose  contour  of  the  hub,  at  some  distance  upstream  of  the  blade 
row,  must  be  replaced  by  some  sort  of  a  smooth  transition).  The  intersection 
of  the  blade  surfaces  with  the  surfaces  n  -  constant  defines  a  two-dimensional 
cascade : 


Figure  2. 


This  cascade  is  then  mapped  into  a  square,  using  the  conformal  transformat, 
described  in  Reference  14.  The  correspondence  of  points  in  the  cascade  a-..: 
mapped  planes  is  shown  in  Fig.  5. 

The  metrics  calculated  in  the  plane  =  constant  must  be  convert  . 
to  the  three-dimensional  quantities  required  by  the  general  coordinate 
transformation,  where  ^  ,  and  L,  are  regarded  as  functions  of  *  ,  r  , 

d  .  If  these  metrics  are  calculated  by  differencing  the  coordinates  t 
selves,  the  easiest  way  to  proceed  is  to  regard  the  sequence  of  plane?  r'  - 
constant  as  determining  f  ,  &  ,  and  X  f°r  given  values  of  %  ,  ^  ,  and  C 
difference  formulas  then  can  use: 

1*  *  (q®*  -  /  £>" 

s*  •  /-a'' 


14.  Rae,  W.J.,  "A  Computer  Program  for  the  Ives  Transformation  in  Turbo¬ 
machinery  Cascades",  Calspan  Report  No.  6275-A-3  (November  IRS^) . 
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where 


*3 


d  C  x  ,  r,  &) 


d($'  *?> 

d(x  ,  r,  0) 

is  the  Jacobian  of  the  transformation. 


rK 

Or 


o. 


(3-20) 


If  the  metrics  in  the  K]  =  constant  plane  are  evaluated  analytically 
(for  example  by  conformal  mapping,  as  in  the  examples  used  here),  then  they 
will  contain  derivatives  of  the  quantity  r&,  taken  at  constant  X  and  n  . 

In  order  to  extract  the  desired  metrics,  the  chain  rule  is  used,  giving: 


££ 

dX 


tl\  +  11) 
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r  90  ). 
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-  11) 

*iJ 


d  n.  /. 
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d  X  /  f 


.*££)♦*)  3AI 

3r/x&  ^  dr'x 


(3- z i ; 


and  similarly  for  the  derivatives  of  £  . 
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Section  4 

BOUNDARY  CONDITIONS 


The  algorithm  described  in  the  previous  section  is  implicit,  in  that 
the  solution  vector  at  all  field  points  is  updated  by  each  sequence  of  three 
one-dimensional  solutions.  The  values  of  the  solution  vector  at  the  grid 
boundaries  can  be  updated  either  as  part  of  this  implicit  scheme,  or  explicitly 
In  the  former  case,  special  finite-difference  versions  of  the  boundary  con¬ 
ditions  must  be  developed,  and  incorporated  into  each  of  the  one-dimensional 
solution  procedures.  In  the  latter,  the  updating  of  the  boundary  values  is 
separated  from  the  sequence  of  finite-difference  operators  used  at  the  field 
points,  and  lags  one  time  step  behind. 

An  explicit  treatment  at  the  grid  boundaries  is  used  in  the  current 
program,  in  order  to  retain  flexibility  with  regard  to  the  coordinate 
transformations  used.  If  an  implicit  treatment  of  the  boundary  conditions 
were  used,  the  specific  details  of  the  transformation  would  have  to  be  built 
into  the  main  solution  algorithm. 


Boundary  conditions  are  needed  for  five  variables:  three  velocity 
components,  and  any  two  thermodynamic  variables,  for  example,  pressure, 
density,  rothalpy,  total  energy).  The  technique  used  in  external-aerodynamic 
studies  (Reference  8,  for  example)  is  to  use  the  surface  tangency  condition  for 
the  three  velocities,  to  extrapolate  the  density  from  nearby  field  points, 
and  to  update  the  pressure  using  an  expression  for  its  normal  derivative  at 
the  surface . 

Sur face-Tangency  Condition 

Since  the  blade  surfaces  lie  in  the  planes  L,  =  const.,  the  surface 
tangency  condition  is 


\ 


£f  -  o  =  it  *  w,  *  w,<r  *  ^  (1-1) 

The  expressions  for  the  other  two  contravariant  components  are  then  added  to 
this,  to  give: 


(f* 
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/ r  \ 

/W*  t 

W,. 

w.  -  <t  \ 

w«  -  ^  j 

Z&/rl 

w*  ) 

-  ?t  J 

Numerical  values  for  W ,  and  are  found  at  the  surface  by  extrapolation, 
after  which  this  equation  is  solved  for  the  surface  values  of  k ,  Wr  ,  and  W-  . 
A  similar  procedure  is  used  at  the  hub  and  tip,  with  W2  =  0  in  that  case. 

Normal  Pressure  -  Derivative  Relation 


The  pressure  at  the  surface  =  constant  is  found  by  using  an 
expression  for  dp/di,  to  extrapolate  from  the  value  a  distance  A  Z,  away 
from  the  surface.  Two-dimensional  counterparts  of  this  derivative  expression 

are  given  in  References  8,  9,  15,  with  relatively  few  details  about  their 
derivation.  The  version  appropriate  to  the  present  problem  is: 


15.  Pulliam,  T.H.  , 
tions  of  Three 


and  Steger ,  J. L. , 
Dimensional  Flow", 


"On  Implicit  Finite 
AIAA  Paper  78-10, 


Difference  Simula¬ 
te  January  1978). 
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1 

I 

This  is  to  be  regarded  as  an  expression  for  d^/dC,  on  the  blade  surfaces 
Z,  =  constant,  with  all  other  quantities  either  known  or  found  by  extrapola¬ 
tion.  Its  usefulness  lies  in  the  fact  that  it  contains  no  time  derivatives 
of  the  dependent  variables.  The  derivation  of  this  equation  is  achieved  by 
summing  the  three  components  of  the  momentum  equation,  multiplied  respectively 
by  ,  £r  and  £8  / r  .  The  resulting  equation  is  then  arranged  into  four 
groups,  containing  derivatives  with  respect  to  z  ,  £  ,7  and  £  respectively. 

Within  each  of  these  groups,  appropriate  terms  are  added  and  subtracted  so 
as  to  form  the  quantity  W3  ,  which  is  zero.  Expansion  of  certain  of  the  ^  , 

7  and  4  derivatives  then  leads  to  cancellation  of  a  number  of  terms  involving 
the  product  of  the  pressure  times  derivatives  of  the  metrics.  Other  terms 
of  this  type,  which  do  not  cancel,  can  be  rearranged  by  noting  the  property 
of  the  Jacobian  that 

1  W  d  3n  d  24 

£>  dec  2Z,  doc  d<[  doc  dt,  doc  1 

where  ct  is  X ,  f  or  9.  After  these  simplifications  have  been  made,  four  of 
the  remaining  eras  can  be  recognized  as  the  continuity  equation;  removing 
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the  dependent  variables.  Finally,  it  is  necessary  to  add  and  subtract  certain 
derivatives  of  /r  and  to  use  the  relation 


r-  d  /±\ ,  d  /J_\  j.  r  3  /J_\  _  3  /  _t_  \  _ 

’s  ^  (rj  '’e  V  r  J  "  (  r> 
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an 
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(4-5) 


A  similar  relation  for  the  pressure  gradient  normal  to  the  hub  and 
shroud  can  be  found,  by  summing  the  momentum  equations,  multiplied  respectively 
by  ,  r\r  ,  and  (and  then  using  f\e  =  0).  The  result  is 


-  -J*  (  'It*  *  uV'l.r) 

+  i  {If  +  ^ 

<■  Pw,  (i*t| 

=  /’Ir  A-®  (4-6) 


At  £  =  -  1,  a  symmetry  condition  is  imposed  ,  explicitly  : 


U  1  .  *[,  $)  "  0  (±  1.  rj  -£) 


(4-") 


Kutta  Condition 


The  trailing-edge  region  is  treated  in  the  present  work  as  though 
the  trailing  edge  were  a  cusp,  i.e.,  the  pressures  and  flow  angles  leaving 
either  side  of  the  trailing  edge  are  required  to  be  the  same,  although  the 
velocity  magnitudes  may  be  unequal.  For  small  trailing-edge  included  angles. 
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the  application  of  the  flow  tangency  condition  on  the  blades  would  be  expected 
to  produce  nearly  equal  flow  angles;  thus  in  the  present  work,  the  pressures 
were  matched,  by  updating  the  quantity  K  at  the  two  points  denoted  by  +  and 
-  in  Fig.  4.  Thus: 

-jo  -  (*-»;  [k  -  ~  (wz-  (uJr)z)~\ 


-P+  -  -  (Tf-  0  [k+-  K~-  ±  ( u*  w*  -  U~  W/*  )] 


=  o 


In  order  to  set  K  and  K  ,  define 


*  *  i 


(4-8) 


(4-9) 


where  the  notation  (  signifies  that  these  are  values  extrapolated  to  the 

surface.  Then  calculate: 


=  (K*-K-)^  -  -L  [U*  wl  -  !j;  U.‘ ]  (4-10) 

-  i<  2  i 
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Figure  4.  Trailing-Edge  Notation 


Conditions  at  Infinity 

The  points  at  infinity  (  ^  =  0,  =  -  1)  are  excluded  from  the  grid 

by  using  an  even  number  of  grid  points  in  the  ^  -  direction: 

Z,=(L-l)AZ>y  1  =  1,2,...,  LMX;  LMX  even,  At,  -  2/LMX-l  (4-11) 

.  +  +  —  — 

The  points  at  ^  =  -  1,  L  =  L  and  L  =  L  (where  L  is  the  integer  part  of 
(LMX  +  10/2,  and  L*  =  L  +  1  -  see  Fig.  3  -  )  are  assigned  fixed  values  at 
the  beginning  of  the  calculation,  and  are  not  changed  thereafter.  In  par¬ 
ticular,  care  is  taken  to  avoid  differencing  across  these  points  when  applying 
the  symmetry  condition  at  =  -  1. 


The  question  of  how  to  select  the  values  that  are  assigned  at  plus 
and  minus  infinity  is  a  serious  problem  in  its°lf.  The  literature  contains  a 


number  of  papers  which  use  time-marching  methods  to  solve  the  Euler  equa¬ 
tions.*6’^  In  several  of  these,  the  method  of  characteristics  is  applied 
to  the  equations  (in  nonconservative  form)  in  order  to  calculate  the  solution 
at  the  grid  boundaries.  In  other  papers,  certain  dependent  variables  such 
as  the  pressure  or  outlet  flow  angle  are  prescribed,  and  the  remaining  vari¬ 
ables  deduced  from  these.  However,  there  does  not  exist  at  present  a  complete 
treatment  of  this  "Trefftz-plane"  problem,  connecting  the  far-field  solution 

to  conditions  at  the  blades,  and  giving  the  relations  between  the  dependent 

->1 

variables  themselves.-  A  prominent  example  of  what  is  missing  from  the 
current  literature .is  the  connection  between  the  pressure  far  downstream  and 
the  trailing-edge  conditions:  in  the  nonlinear  small-disturbance  theory, 
the  imposition  of  the  Kutta  condition  at  the  trailing  edge  uniquely  deter¬ 
mines  the  circulation  at  that  spanwise  station,  which  in  turn  defines  the 

pressure  rise  and  turning  angle  that  must  be  reached  far  downstream  of  the 
1  4 

blade  row.  ’  The  extension  of  this  relationship  to  the  case  of  the  full 
Euler  equations  has  not  been  made.  Thus,  for  example,  it  must  be  presumed 


16.  McDonald,  P.W.,  "The  Computation  of  Transonic  Flow  Through  Two-Dimen¬ 

sional  (las  Turbine  Cascades",  American  Societv  of  Mechanical  Engineers, 
Paper  71-GT-89,  1971. 

17.  Gopalakrishnan,  S.  and  Bozzola,  R. ,  "Computation  of  Shocked  Flows  in 

Compressor  Cascades",  .American  Society  of  Mechanical  Engineers,  Paper 
''2-CT-31,  1972. 

18.  Kurzrock,  J.W.  and  Novick,  A.S.,  "Transonic  Flow  Around  Rotor  Blade 
Elements",  Transactions  of  the  .American  Society  of  Mechanical  Engineers, 
Vol.  97,  December  197S,  pp.  598-607. 

19.  Thompkins,  William  T. ,  Jr.,  ".An  Experimental  and  Computational  Study 
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Power  99,  (1977)  53-62. 

21.  Karamcheti,  K. ,  Principles  of  Ideal-Fluid  Aerodynamics,  Wiley  and  Sons, 
New  York  (1966)  Section  19.4. 
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that  the  assignment  of  an  outlet  flow  angle  implies  a  violation  of  the  Kutta 
condition,  and  this  in  turn  renders  the  circulation  and  pressure  rise  non¬ 
unique,  It  is  not  clear  from  the  published  literature  how  these  problems  are 
resolved  in  current  computer  codes. 

The  present  research  has  not  addressed  these  questions,  because  of 
the  instabilities  encountered  in  the  process  of  developing  the  program.  In¬ 
stead,  a  set  of  boundary  values,  described  below,  was  assigned  at  downstream 
infinity.  These  values  were  adequate  for  use  in  the  early  development  of  the 
computer  program,  but  will  have  to  be  replaced  by  a  more  exact  formulation, 
after  the  grid-related  oscillations  have  been  removed. 

The  specific  choices  for  the  variables  at  downstream  infinity  were 
made  as  follows:  the  static  pressure  ratio  across  the  blade  row  was  assigned 
as  an  input,  and  the  density  was  found  from  the  isentropic  relation.  The 
area  ratio  between  outlet  and  inlet  was  assigned,  and  from  this  the  axial 
velocity  component  was  found  by  conserving  the  mass  flow  pR The  radial 
velocity  Wr was  set  equal  to  zero,  and  the  circumferential  component  W9  was 
chosen,  following  Reference  19,  as  that  value  which  gives  a  uniform  static 
pressure,  i.e.,  radial  equilibrium  requires 

f  ^7  =  °  =  "r  ^9  ~UJrlZ>  or 

The  value  of  the  quantity  K  then  follows  from  these  specifications. 

Wake  Conditions 


Whenever  a  compressor  blade  is  acted  on  by  a  lift  that  varies  with 
radius,  a  sheet  of  vorticity  will  be  shed  from  the  trailing  edge.  The  origins 
of  this  vortex  sheet  can  be  seen  in  the  nonlinear  small-disturbance  results 
of  Reference  4;  Figure  11  of  that  paper  shows  the  distributions  of  radial 
velocity  at  a  sixty-percent  chord  location.  These  distributions  retain  the 
same  qualitative  behavior  all  the  way  to  the  trailing  edge,  i.e.,  they  reveal 
a  discontinuity  at  the  trailing  edge,  which  trails  downstream  as  a  vortex 
sheet . 
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In  the  small-disturbance  potential-theory  approximation,  the  trailing 

vortex  sheets  are  assumed  to  lie  on  the  helical  surfaces  defined  by  the  inlet 

flow.  In  a  full  nonlinear  treatment,  they  must  be  allowed  to  deform,  away  from 

these  surfaces,  as  they  move  downstream.  This  problem  has  been  studied  re- 

22-24 

cently  in  a  series  of  papers  by  McCune  and  Hawthorne.  These  papers 

constitute  a  basis  on  which  to  model  the  vortex-sheet  trajectories  in  a 
finite-difference  code.  For  example,  the  discontinuities  in  radial  velocity 
that  occur  at  the  trailing  edge  would  have  to  be  inserted  at  the  trailing-edge 
location,  and  the  locus  of  this  discontinuity  would  have  to  be  followed  down¬ 
stream.  The  conformal  transformation  used  in  the  present  work  is  not  well 
suited  for  doing  so,  however,  since  the  path  followed  by  the  trailing  vortex 
sheet  in  the  %  ,  4  plane  (Figure  3)  is  a  line  that  leaves  the  trailing-edge 
image,  arid  spirals  around  the  image  of  the  point  at  downstream  infinity.  It 
would  be  virtually  impossible  to  use  a  grid  that  is  fine  enough  to  resolve 
these  discontinuities  numerically.  In  order  to  facilitate  the  resolution  of 
the  vortex-sheet  behavior,  it  would  be  necessary  to  use  a  different  coordinate 
transformation,  in  which  the  path  of  the  vortex  sheet  is  not  as  convoluted 
as  it  is  in  the  case  of  the  Ives  mapping. 

As  an  alternative,  it  might  be  possible  to  trace  the  magnitude  and 
location  of  these  discontinuities  by  "floating  vortex-sheet  fitting",  in 
analogy  to  the  procedure  of  floating  shock  fitting.  The  development  of  such 
a  procedure  was  not  considered  in  the  present  research. 

22.  McCune,  J.E.,  and  Hawthorne,  W.R.,  "The  Effects  of  Trailing  Vorticity 
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Section  5 
RESULTS 


A  two-dimensional  cascade  was  chosen  for  the  purpose  of  debugging 
the  program.  The  blade  geometry,  show™  in  Figure  5,  is  the  same  as  used 
by  Rae  and  Homier:  (Reference  25).  It  was  not  chosen  on  the  basis  of  any 
design  method,  but  only  for  the  purpose  of  facilitating  a  demonstration  cal¬ 
culation.  The  solidity  is  moderate,  and  the  large  leading-edge  radius 
minimises  strong  flow  field  gradients  in  that  region. 


In  order  to  do  a  two-dimensional  case,  the  hub  and  shroud  radii  were 
taken  to  be  t*  /Ca  =  99.5,  •  Cr  / =  100.5,  and  the  calculations  were  done 
at  f / Cg^  =  100.  The  inlet  relative  Mach  number  and  flow  angle  were  taken 
as  0.5  and  33°.  All  quantities  were  made  dimensionless  by  dividing  by  the 
appropriate  combination  of  the  density  and  axial  velocity  component  far 
upstream,  w  and  U_oa  ,  and  the  axial  projection  of  the  chord,  .  Thus, 
for  example,  the  dimensionless  angular  velocity  was  input  as 


^C-q,  __  wr  .  _  O'-oo 

U  -  oa  O-.co  r  ^-oo 


A. 


i 

TOO 


? 


TOO 


(5-1) 


The  specific-heat  ratio  was  taken  as  1.4,  and  the  grid  sices  in  the  transformed 
plane  as  /<MX=  11,  LMX  =  10,  giving  =  0.2,  AZ,  =  2/9. 


To  start  the  calculations,  an  input  tape  was  prepared,  containing 
values  of  the  metrics,  the  Jacobian,  and  the  radii  at  each  grid  point.  On 
the  first  run,  all  dependent  variables  were  initialized  at  their  upstream 
values.  At  the  end  of  each  run,  the  metrics,  Jacobians,  and  radii  were 


25.  Rae,  W.J.,  and  Homicz,  G.F.,  "A  Rectangular-Coordinate  Method  For  Cal¬ 
culating  Nonlinear  Transonic  Potential  Flow-fields  in  Compressor  Cascades", 
AIAA  Paper  78-248  (January  1978). 
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rewritten  on  a  new  tape,  along  with  values  of  the  solution  at  the  last  tmf 
step.  This  tape  could  then  be  used  to  start  the  next  series  of  time  steps. 
The  maximum  time  step  allowed  by  the  CFL  condition  (see  Section  3)  was 
calculated  at  the  end  of  each  step;  all  of  the  results  discussed  below  we  re¬ 
calculated  with  a  time  step  equal  to  half  this  value,  and  with  damping  coef 
ficients  and  Eg  equal  to  the  (dimensionless)  value  of  the  time  step.  Th 
static  pressure  ratio  was  assigned  as  1.1,  and  the  area  ratio  as  1.0. 

Figure  6  shows  results  for  the  velocity  field  after  ten  time  steps 
In  the  guided  channel  between  the  blades,  the  results  conform  to  what  w ■•■.Id 
be  expected,  but  at  the  stations  K  =  2  and  K  =  10,  there  are  very  large  os-; 
lations.  On  this  coarse  grid,  these  two  stations  are  the  outermost  ones  a: 
which  implicit  calculations  are  done;  the  stations  K-  1  and  11  form  the 
boundaries  of  the  computational  grid,  and  their  values  are  updated  exp  lie:* 
one  time  step  behind. 

The  oscillations  at  the  station  K  =  2  are  due  to  the  metric  singu¬ 
larities  at  K  =  1;  the  image  of  the  point  at  upstream  infinity  is  at  X  = 
and  midway  between  the  central  pair  of  L  -  values  (see  Figure  3).  Even 
though  the  metrics  at  K  -  1,  L  -  L+  and  L  are  finite,  nevertheless  their 
gradients  are  so  steep  at  those  points  that  they  destabilize  the  solution. 
This  effect  can  be  displayed  clearly  by  the  results  of  a  calculation  in  w... 
the  solution  is  initialized  to  the  freestream  values,  and  then  advanced  - 
a  single  time  step.  For  the  two-dimensional  case,  the  basic  equation 
becomes 
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But  the  quantities  in  the  square  brackets  on  the  right  side  of  this  equation 
are  each  zero  (as  can  be  verified  from  Eq.  3-20,  with  ?r/ dr\  =  1).  This  result 
is  exact,  analytically,  but  when  evaluated  numerically,  especially  on  a  coarse 
grid,  the  result  is  nonzero,  and  of  a  magnitude  consistent  with  the  magnitude 
of  the  oscillations  that  develop.  Part  of  the  problem,  in  the  present  case, 
is  due  to  the  use  of  analytic  formulas  for  the  metrics;  it  is  pointed  out  in 
Reference  8  that  metrics  which  are  generated  numerically,  by  differencing 
the  coordinate  mapping  itself,  are  actually  less  sensitive  to  this  problem 
than  the  analytic  metrics. 

In  an  attempt  to  alleviate  this  problem,  the  values  of  the  metrics 

at  K-  1,  L  -  L'  and  L+  were  changed,  as  follows:  the  value  of  £„/r) 

L  -l+  9  * ' 1 

which  is  equal  to  -  <£a/r)  ^  ■>  was  chosen  so  as  to  make  the  factor  multiply- 

ing  G<o  numerically  equal  to  zero.  Mext ,  the  value  of  £  J  *f  ,  which  is 

L  =  L.* 

equal  to  -  £ )  ,  was  changed  so  that  the  factor  multiplying  £  would 

have  the  same  value  at  K  =  1,  L  -  L~  and  L+  .  The  result  of  this  modifica¬ 
tion  is  shown  in  Figure  7,  at  the  tenth  time  step.  Comparison  with  Figure  6 
shows  that  a  considerable  smoothing  of  the  flow  pattern  at  K  =  1  was  achieved. 

Analogous  modifications  were  made  at  K  =  KMX,  but  these  did  not 
yield  the  same  degree  of  success,  presumably  because  the  solution  on  this 
line  is  also  affected  by  the  singular  region  near  the  trailing  edge.  Several 
attempts  were  made  to  overcome  this  problem,  by  using  extrapolation  to  update 
the  points  near  the  trailing  edge.  These  attempts  were  not  successful.  More¬ 
over,  this  approach  is  difficult  to  justify,  since  points  on  the  surface 
near  the  trailing  edge  are  very  important  to  the  solution,  in  that  they  are 
the  ones  used  in  applying  the  Kutta  condition,  as  well  as  in  enforcing  the 
surface  tangency  condition.  Any  extrapolation  procedure  alters  the  role  of 
these  boundary  points,  making  them  dependent  on  the  field  behavior,  rather 
than  the  other  way  around. 

The  metric-singularity  problems  encountered  in  this  research  are 
aggravated  by  the  use  of  the  coarse  grid,  and  by  the  fact  that  the  grid  wraps 
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around  the  trailing  edge  of  the  blades.  While  the  use  of  a  finer  grid  might 
alleviate  the  situation  somewhat,  it  was  not  used,  since  it  appeared  that 
lengthier  calculations  would  not  be  justified  until  the  more  basic  cause  - 
the  trailing-edge  singularity  -  was  eliminated. 
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CONCLUDING  REMARKS 


The  experience  gained  in  this  study  suggests  strongly  that  type  of 
grid  used  (where  the  metrics  are  evaluated  from  analytic  formulas,  and  the 
grid  is  wrapped  around  the  trailing  edge)  is  not  suitable  for  use  with  the 
implicit  time-marching  algorithm.  The  evidence  is  not  conclusive,  however; 
it  may  be  that  use  of  a  finer  grid  and  some  other  special  treatment  of  the 
metrics  in  the  region  of  the  singularities  could  stabilize  the  calculations. 
However,  a  preferable  course,  for  future  developments,  appears  to  be  the  use 
of  grids  which  are  free  from  singularities,  especially  in  the  trailing-edge 
region. 


After  the  grid-induced  instabilities  have  been  removed,  a  number  of 
other  problems  will  remain.  These  problems  were  not  considered  in  depth 
during  this  research,  because  of  the  amount  of  effort  devoted  to  the  insta¬ 
bility  problems.  .Among  the  topics  that  will  have  to  be  considered  are  the 
Kutta  and  far-field  conditions,  the  location  and  strength  of  the  trailing  vor¬ 
tex  sheets,  and  the  interaction  between  all  of  these.  Also  remaining  is  the 
problem  of  shock  capturing  in  genuinely  transonic  flows,  which  may  require 
alterations  in  the  difference  formulas  used.  Finally,  there  are  several  ad¬ 
vances  in  the  time-marching  algorithm  that  have  taken  place  during  the  period 
of  this  research,  such  as  the  technique  of  flux  vector  splitting. “  These 
should  be  considered  for  incorporation  in  the  numerical  method. 


26.  Steger,  J.L.,  and  Warming,  R.F.,  "Flux  Vector  Splitting  of  the  Inviscid 
Gasdynamic  Equations  with  Application  to  Finite  Difference  Methods", 

NASA  TM  7S605  (July  1979). 
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The  dependent  variables  are  defined  as: 
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Thus  the  matrices  A,  B,  C,  and  D  are 
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